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Simulating the Common Envelope Phase of a Red Giant Using 

SPH and Uniform Grid Codes 

Jean-Claude Passy 1,2 , Orsola De Marco 3 , Chris L. Fryer 4 , Falk Herwig 2 , Steven Diehl 4 , 
Jeffrey S. Oishi 5 , Mordecai-Mark Mac Low 1 , Greg L. Bryan 6 , Gabriel Rockefeller 4 

ABSTRACT 

We use three-dimensional hydrodynamical simulations to study the rapid in- 



fall phase of the common envelope interaction of a red giant branch star of mass 
equal to 0.88 M & and a companion star of mass ranging from 0.9 down to 0.1 M . 

^vq We first compare the results obtained using two different numerical techniques 

with different resolutions, and find overall very good agreement. We then com- 

P^ pare the outcomes of those simulations with observed systems thought to have 

gone through a common envelope. The simulations fail to reproduce those sys- 
tems in the sense that most of the envelope of the donor remains bound at the end 
of the simulations and the final orbital separations between the donor's remnant 
and the companion, ranging from 26.8 down to 5.9 R & , are larger than the ones 
observed. We suggest that this discrepancy vouches for recombination playing 
an essential role in the ejection of the envelope and/or significant shrinkage of 

CN the orbit happening in the subsequent phase. 

cs 

Subject headings: binaries: close — binaries: general — hydrodynamics — meth- 
ods: numerical — stars: evolution 



1. Introduction 



Around 60% of F and G stars are binaries, of which about 30% have separations smaller 



than 30 AU and will interact during the primary's evolution (Duquennoy & Mayor 1991). 
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During the giant phases of the primary, companions closer than ~ 5 AU enter a strong 
interaction phase with the primary and, under certain circumstances, a common envelope 
(CE) may form around the two stars. The secondary star spirals inside the envelope of 
the primary and may also fill its own Roche lobe because it cannot accrete all the matter 
coming from the donor star. This process is called a common envelope interaction and was 



originally described by Paczynski (1976). For a general review of the topic see, e.g., Iben 



& Livio (1993). There are two different processes leading to the onset of a CE phase: the 



start of unstable mass transfer from the expanding primary to the secondary (Hjellming & 



Webbink 1987; Hurley et al. 2002) and the development of a tidal instability that occurs if 



there is not enough angular momentum in the orbit to maintain the primary's envelope in 
synchronization ( Darwin|1879 ). The post-CE system will be either a compact binary system, 
if there is enough energy to eject the primary's envelope, or a merger, if not. 

The CE interaction is an essential ingredient for any binary population synthesis study of 
intermediate (e.g., Politano et al.poTO ) or massive stars (e.g., Belczynski et al.||2008 ). Com- 
pact binaries are believed to be formed through at least one CE phase. Among them are 
symbiotic binaries, supersoft X-ray sources, cataclysmic variables and double white dwarfs, 



which are all possible supernova Type la progenitors. As Meng et al. (2010) pointed out 



results deduced from population synthesis studies such as the Type la supernova birth rate 
are highly dependent on the physics of the CE phase. Therefore, it is paramount to under- 
stand more accurately the CE interaction in order to identify the formation channels of such 
supernovae and to compare observations with predictive models. Moreover, many substellar 
companions to evolved stars have recently been discovered with small orbital separation. 



Maxted et al. (2006) found a brown dwarf orbiting a white dwarf with a 116 min period, 



while Setiawan et al. (2010) discovered a system composed of a Jupiter-like object orbiting 



an horizontal branch star with a 16.2 days period. We therefore know that substellar com- 
panions can survive a CE interaction, but what is the minimum mass of the companion that 
can eject the envelope? Is the ejected envelope entirely unbound or will some of it eventually 
fall back and form a circumbinary disk? Were the substellar companions present before and 
survived the CE or were they formed later on in such a disk ( PeretspOlO )? Those questions 
remain unanswered. 

Although the CE process was outlined more than 30 years ago, it is still far from 
understood quantitatively. Numerical simulations suggest that the typical duration of the 
entire CE phase is short - - less than 10 3 years - - which makes CE ejections unlikely to 
be observed. However, one can use observations of post-CE binaries to better understand 
CE evolution. With the use of stellar models, the initial configuration of such systems can 
be approximately determined from the final configuration. Using either the a-formalism 



(Webbink 1984, but see De Marco et al. 2011) or the 7-formalism (Nelemans et al. 2000) 



the relevant parameters can be constrained and the CE ejection efficiency can be predicted. 



Using this approach, De Marco et al. (2011) suggested an anti-correlation between a, the 



CE efficiency parameter, and the secondary to primary mass ratio. 

The entire CE evolution can be divided into three different phases ( Podsiadlowski|2001 ) 
with different timescales, length scales and physics involved. These differences are the rea- 
sons why reproducing the entire CE evolution of a given system accurately is challenging. 
Therefore, one usually treats one phase after the other with different methods. In this paper, 
we focus only on the rapid infall phase, which has a short timescale (~ 1 — 10 years), and 
in which the evolution is driven by drag forces. Several numerical hydrodynamic studies of 



the CE interaction have been carried out in the past (for an exhaustive list, see Taam & 



Sandquist 2000), including a series of ten papers starting with the two-dimensional calcu- 



lation of the interaction of a 16 M supergiant and a 1 M neutron star (Bodenheimer & 



Taam|[l984 ) , and most recently treating three-dimensional simulations of the CE interaction 



between 3 or 5 M giant stars and 0.4 or 0.6 M main sequence (MS) companions (Sandquist 
et al.|1998 ). The latter study has been extended first by Sandquist et al. (2000) to 1 M and 
2 M Q red giant branch (RGB) stars with companion masses ranging from 0.1 to 0.45 M , 



then by De Marco et al. (2003) to a 1 M asymptotic giant branch star with a 0.1 or 0.2 M 6 



© 



companion. Ricker & Taam (2008) computed high resolution simulations of the CE phase 



between a 1.05 M RGB star and a 0.6 M compact companion, and concluded that the 
gravitational component of the drag dominates over the hydrodynamical component (also 

see 



Taam & Ricker 2010 Ricker & Taam 2011). 



A direct comparison of the results obtained using different numerical methods has how- 
ever never been carried out. Although analytical/empirical work has included discussion 
regarding observational data, there are only a couple of publications that connect simula- 



tions and observations in a meaningful way (see e.g., Sandquist et al. 2000). Those are, as 
we will explain in 



J and £J4J key steps to better understand the implications of CE interac- 
tions and the physical processes driving them. In this paper we therefore present numerical 
simulations with two different algorithms of the CE interaction of a 0.88 M RGB star with 
a MS companion. Different companion masses from 0.1 M to 0.9 M are considered. The 



simulations are carried out with both an Eulerian code (Enzo in uniform-grid mode, O'Shea 



et al. (2004) and enzo.googlecode.com) and a Lagrangian code (SNSPH, Fryer et al. 2006), 
and for different resolutions. We describe the numerical methods and the initial conditions 
of our 15 simulations in £[2] and §[3j We describe and discuss the results in ^4] and §[5j and 
finally conclude and summarize in $6} 



2. Description of codes 

In this section we describe the numerical methods we use. We first compare the code 
algorithms and explain why a code-to-code comparison is necessary. Then, we describe both 
codes in detail and finally discuss different ways to compare resolution. 



2.1. Eulerian vs Lagrangian codes 

Although they are meant to simulate similar astrophysical situations, high order Eule- 
rian grid codes and Lagrangian smoothed-particle hydrodynamics (SPH) codes differ fun- 



damentally, with each having advantages and disadvantages. Among other studies, Davies 



et al. (1993), Frenk et al. (1999), Agertz et al. (2007), Tasker et al. (2008) and Heitsch 



et al. (2011) aim at identifying these differences. On the one hand, high-order Eulerian grid 



codes have a better wavenumber resolution than SPH codes for an equal number of cells and 
particles and are more accurate at resolving the rarefied regions since, unlike SPH, the res- 
olution does not depend on the density of the gas; Eulerian codes also better resolve shocks 
( Tasker et al.|2008 ) compared to SPH codes; and finally, SPH noise dominates subsonic flows 
and therefore makes it difficult for SPH codes to follow perturbations in flows with Mach 
numbers under unity. On the other hand, SPH codes don't diffuse material properties, and 



inherently conserve mass, momentum and energy (Rosswog 2009). While the treatment of 



boundary conditions can be challenging in grid-based codes when the flow expands beyond 
the computational domain, SPH easily handles vacuum conditions. It is still unclear which 
method is the most appropriate to simulate CE interactions. Therefore, we use both methods 
and confront the results from both codes in order to draw conclusions about their physical 
relevance. 



2.2. Input physics 

Both codes solve the fully compressible hydrodynamics equations with self-gravity in- 
cluded. These equations can be written using an Eulerian formulation: 
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where p, v,p, $,M,M ini ,7 are the density, velocity, pressure, gravitational potential, specific 
total energy, specific internal energy and adiabatic index of the gas, respectively. The total 
energy is the sum of the internal energy, and the macroscopic kinetic energy: 



u 



U;, 



mt 



+ v 2 /2. 



(6) 



Equations (HI), $2h and pi) express mass continuity, conservation of momentum and conser- 
vation of energy, respectively. Both codes evolve the internal energy rather than the total 
energy. An ideal gas equation of state (Eq. [I]) for a monoatomic gas (7 = 5/3) closes the 
system composed by equations Q-(|3]). Such an equation of state represents an adequate 



approximation of the deep convective envelope of RGB stars (Hjellming & Webbink 1987) 



although it ignores some physical processes such as radiation pressure and ionization. We 
discuss this point in detail in §5.2.2 Finally, the gravitational potential is calculated using 
the Poisson equation (Eq. [5]). 



2.3. The Enzo code 
Enzo is a three-dimensional, adaptive mesh refinement hybrid (hydrodynamics + N- 



body) grid-based code (Bryan et al. 1995; O'Shea et al. 2004) that we use in uniform-grid 



mode only. It is primarily designed to simulate cosmo logical structure formation (Norman 



et al. 2007). However, its numerous features make it useful for reproducing many different 



astrophysical situations, including CE interactions. 



The Euler equations (Eqs.nTJ3l are solved using the van Leer (1977) second-order advec- 



tion method also implemented in Zeus (Stone & Norman 1992). Although those equations 



can also be solved in Enzo by a third-order piecewise parabolic method that better resolves 
shocks and turbulence, our tests show that it slows down the computation by a factor 2. 
As we will point out in $4l there are neither strong shocks nor important turbulence in our 



simulations so we favor efficiency and use the van Leer solver. The Poisson equation is solved 
using fast-Fourier transforms. 

In the case of a CE interaction between a RGB star and a MS companion, the radius 
of the secondary -- typically 0.5 R & --is small compared to the primary's radius (~ 100 
Rq), so we can legitimately model the companion as a point mass particle. Furthermore, as 
shown in Fig. [31 the primary's core is also small (~ 0.01 R Q ) and dense, so it can also be 
modeled as a point mass. 

Enzo usually models collisionless particles as a continuous mass field appropriate for 
computing the gravitational potential in the case that each particle represents many actual 
particles, such as in cosmological simulations with dark matter. In that case their mass is 
deposited in the 8 nearest cells and added to the gas density of those cells to find the total 
density for use in solving the Poisson equation (Eq. [5]). In a simple two-body interaction 
between 1 M and 0.1 M© objects in a one year circular orbit without gas, this method does 
not provide the accuracy required by our problem because of the spreading out of the mass 
of the point source, leading to an inaccurate gravitational potential. Indeed, a 1 % error in 
the orbit is reached after only 6 orbits. Consequently, we implemented, as a new type of 
particle, point mass particles. These particles create a potential that is added analytically to 
the gas potential calculated using the Poisson equation. Using an analytic potential yields 
an accuracy of the orbit more than two orders of magnitude better than with the default 
particles. The gravitational potential created by a point mass particle is smoothed according 



to the prescription of Ruffert (1993), used in Sandquist et al. (1998): 



^ i \ -GM PM 

yr 2 + e 2 5 2 exp [— r 2 /(e<5) 2 ] 

where Mpm is the mass of the particle, r is the distance from the particle, 6 is the size of a 
cell and e = 1.5. The point mass particles are advanced using a leapfrog algorithm. Time 
stepping is determined by taking the minimum time step between the Courant conditions 
for the gas, the particles and the acceleration field: 
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where C\ = 0.4 is the Courant factor, C2 = 0.4 is the particle Courant factor, c s is the sound 



speed, v 



[V X , Vy, V Z/ 



is the velocity of the gas, V = (V x , V y , V z ) is the velocity of a particle 



and g = (g xi g y ,g z ) is the acceleration field. 

Finally, we remark that the current Enzo Poisson solver prevented us from using nested 
or adaptive grids that would have allowed us to increase resolution locally. The inaccurate 
treatment of boundary conditions within the refined grids prevented us from stabilizing the 
RGB progenitor in a multi-grid initial setup. We are currently developing a new Poisson 
solver that will allow us to use nested grids as well as adaptive mesh refinement and carry 
out better-resolved simulations. 



2.4. The SNSPH code 



SNSPH (Fryer, Rockefeller, & Warren 2006) is a three-dimensional, parallel SPH code 



using tree gravity. It uses a regular Monaghan cubic spline kernel (Monaghan 1992). For 



the artificial viscosity we use the sum of a bulk viscosity and a von Neumann and Richtmyer 



viscosity (Rosswog 2009). The particles are organized into a parallel hashed oct-tree as 
described in Warren & Salmon (1993). The gravitational potential of a SPH particle, i, is 



smoothed using the following formula: 



®i(Xi = Ti/ hi) 



( -Grm/hi x (|x? - ^xf + ±xl - 1.4) if < x < 1 

-Grm/n x [(|x t ? - xf + ±xf - is? - 1.6) /hi + l/15rj if 1 < x < 2 



-Grrii/ri 



otherwise 
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where hi, rrti and r, are the smoothing length, the particle mass and the distance from the 
particle, respectively. We compare both numerical potentials to the theoretical potential in 
Fig. [11 For a given smoothing length, hi, the Monaghan (1992) potential used in our SNSPH 



simulations is deeper than the Ruffert (1993) one used in our Enzo simulations. Also, the 



Monaghan (1992) potential is exact at distances larger than 2h{ whereas the Ruffert (1993) 



potential only asymptotically tends to the exact potential. 



SNSPH uses the fast multipole method to calculate gravitational accelerations (Warren 



fc Salmo"nfl i~993 ) . The SPH particles are also advanced using an leapfrog algorithm. Finally, 
in order to keep the same overall spatial coverage, the smoothing length varies according to 



the formula from Benz (1989): 
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Fig. 1. — Comparison between the different potentials in arbitrary units with hi = eS = 1. 
Plotted are the theoretical potential (solid line), the Ruffert ( |1993 ) potential used in Enzo 



(dashed line) and the Monaghan ( 1992 ) one used in SNSPH (dash-cross line) 



2.5. Resolution comparison 

There is no ideal way to compare the resolution between SPH and uniform-grid codes. 
However, a few criteria can give us a general idea of how to relate them. 



As mentioned by Davies et al. (1993), a first global criterion would be to compare the 



total number of SPH particles N part with the total number of cells originally inside the 
progenitor: 
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where N tot , V\, Vq, R\ and L are the total number of cells, the volume of the primary, 
the volume of the grid, the radius of the primary and the linear dimension of the grid, 
respectively. As time goes by, the gas will however fill a larger fraction of the numerical grid 
and thus increase the number of relevant cells, but not the real resolution of the simulation. 

A more local criterion is to compare the size of an Enzo grid cell, 5, with the SPH 
smoothing length, which varies in space and time. Indeed, if the companion does not sink 
much into the primary's envelope and does not modify the inner part of the smoothing 
length distribution too much, then the resolution deep inside the progenitor does not matter. 
Therefore, we compare the smoothing length distribution of the SPH model to the cell size 
of the Eulerian grid. As shown in Fig. |2j the smoothing length at small radii does not vary, 
so an Enzo run with a 128 3 grid will be under-resolved compared to our canonical 500 000 
(roughly 80 3 ) particle SPH run no matter how deep the companion penetrates while a run 
with a 256 3 grid would be equivalent to our SPH runs if the separation between the primary 
core and the companion always exceeds 20 R & . This local criterion for the resolution is 
not perfect either since it does not take into account the variation of the smoothing length 
throughout the SPH simulation. 
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Fig. 2. — Resolution comparison between the SNSPH smoothing length field (dots) for a run 
with 500 000 particles, and the Enzo size of a grid cell for the 128 3 (dash line) and the 256 3 
(solid line) runs for the initial (left) and final (right) particle distributions. 
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Again, comparing the resolution between uniform-grid and SPH codes is quite challeng- 
ing and both methods have, in the situation we are interested in, strengths and weaknesses: 
SPH will under-resolve the low-density outer parts of the envelope, where the smoothing 
length dramatically increases, while it will be more accurate in the later phase of the evolu- 
tion when the separation between the primary's core and the secondary will typically sink 
below a few cells. Therefore, comparing SPH and grid-based simulations is paramount in 
order to state which one is more adapted to our problem, and the combination of both global 
and local criteria is the best way to compare the resolutions of both methods. 



3. The simulations 

We perform 5 SNSPH and 12 Enzo simulations of CE interactions with a 0.88 M & 
RGB primary that are summarized in Table [TJ The SNSPH simulations are computed using 
500 000 particles whose initial smoothing length follows the radial profile shown in Fig. [2] 
The Enzo simulations are performed using either a 128 3 or a 256 3 grid. In both cases, the 
linear size of the computational domain is L = 3 x 10 13 cm. We consider companion masses 
of 0.9, 0.6, 0.3, 0.15 and 0.1 M . Giant stars are slow rotators with rotational velocities of 



the order of a few kms 1 (de Medeiros & Mayor 1999). Although it is expected that a close 



companion will, through the action of tides and the transport of angular momentum in the 
primary envelope, spin up the envelope during the pre-CE phase, the actual rotation of the 
primary at the onset of the CE interaction is hard to quantify. Moreover, even if the primary 
was uniformly rotating at 50 kms -1 , its rotational energy would be 

E rot = -rgMiRJcu 2 ~ 2.2 x 10 44 ergs (14) 

where u, r g , Mi and Ri are the angular velocity, the radius of gyration, the mass and 



the radius of the primary, respectively. For RGB stars r g is typically about 0.1 (Taam & 
Sandquis~t||2000 ) . This rotational energy does not affect the energetics of the system since it 
is more than two orders of magnitude smaller than the binding energy of the primary (see 
below). Consequently, we assume that the primary is initially non-rotating. Finally, the 
companion is at the start placed at the surface of the primary in a circular orbit. We thus 
have three different simulations for each initial companion mass - one with SNSPH, and two 
with Enzo on 128 3 and 256 3 grids. Additionally, we also run two Enzo simulations in order 
to study the dependency of the final parameters on the initial conditions. We consider the 
128 3 Enzo simulation with a 0.3 M companion (Enzo3) as the reference and run identical 
simulations increasing, by 5 %, either the initial velocity of the companion (Enzoll) or 
the initial separation (Enzol2). All the runs follow the evolution of the system for about 
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1 000 days. 





Npart Or NoelU 


M 2 (M ) 


A (R Q ) 


P (days) 


vo/v circ 


A f (R e ) 


P f (days) 


SPH1 


500 000 


0.9 


83 


66 


1 


26.8 


13.5 


SPH2 


500 000 


0.6 


83 


72 


1 


20.6 


10.1 


SPH3 


500 000 


0.3 


83 


81 


1 


11.3 


5.5 


SPH4 


500 000 


0.15 


83 


86 


1 


7.3 


3.0 


SPH5 


500 000 


0.1 


83 


88 


1 


6.1 


2.2 


Enzol 


128 3 


0.9 


91 


75 


1 


28.1 


15.5 


Enzo2 


128 3 


0.6 


91 


83 


1 


20.0 


11.0 


Enzo3 


128 3 


0.3 


91 


93 


1 


11.7 


5.6 


Enzo4 


128 3 


0.15 


91 


99 


1 


8.6 


3.4 


Enzo5 


128 3 


0.1 


91 


102 


1 


8.5 


3.3 


Enzo6 


256 3 


0.9 


85 


68 


1 


25.5 


13.2 


Enzo7 


256 3 


0.6 


85 


75 


1 


19.2 


9.8 


Enzo8 


256 3 


0.3 


85 


84 


1 


11.2 


5.4 


Enzo9 


256 3 


0.15 


85 


89 


1 


6.9 


2.8 


EnzolO 


256 3 


0.1 


85 


92 


1 


5.7 


2.1 


Enzol 1 


128 3 


0.3 


91 


93 


1.05 


12.0 


4.6 


Enzol2 


128 3 


0.3 


95.5 


99 


1 


12.2 


5.0 



Table 1: Main parameters for the different simulations. Reported are the number of particles 
{Np ar t) or cells (Nceiis), the companion mass, the initial orbital separation (Aq), the initial 
orbital period (Pq), the ratio of the initial orbital velocity of the companion (vq) to the 
velocity required for a circular orbit (v circ ) and the final orbital separation (Af) taken at the 



end of the rapid infall phase (§ 4.1). 



As a primary, we use a one-dimensional model of a star with a MS mass of 1 M . Using 



the stellar evolution code EVOL (Herwig 2000), this progenitor was evolved to the RGB 
phase until the core reached M c = 0.392 M . At that time, the radius of the star was 83 R Q 
and its total mass was M 1 = 0.88 M due to mass loss, which was treated using the Reimers 
formalism with r) = 0.5. We adapt this model by using the density and pressure profiles, 
but computing the internal energy using Eq. |4j A sample of relevant profiles are plotted in 
Fig-i 

We now explain how this stellar model is modified in order to be compatible with 
an input suitable for each of our codes. For the SNSPH simulations, the initial particle 



configuration is a weighted Voronoi tesselation (WVT) similar to that described by Diehl & 



Statler (2006). As we have explained in £2.1 one limitation of SPH codes is the large number 



of particles required by dense regions such as the core of the primary. Since the time step 
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Fig. 3. — Comparison between the EVOL stellar evolution model (blue), the SPH initial 
model computed with 500 000 particles (black) and the Enzo initial models for a 128 3 (green) 
and a 256 3 (purple) unigrid. The vertical line represents the core-envelope boundary accord- 



ing to the criterion of De Marco et al. (2011). 



induced by a particle i can be roughly estimated by hi/c S} i where c S} i is the local sound speed, 
a small smoothing length will require a small time step resulting in a high computational cost. 
Since the equation of state changes significantly around the helium core, we represent the 
core by a particle with mass M c . The associated smoothing length is h c = 0.1 R & . We add 
SPH particles in the region around the core such that the density values and gradient profiles 
connect smoothly at the core/envelope boundary (r = 2h c = 0.2 R®), In this way, we obtain 
the profile shown in Fig. [3j Since the density profile has been changed, one must modify 
the gravitational acceleration accordingly. Assuming hydrostatic equilibrium in spherical 
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symmetry, we integrate the pressure gradient choosing the integration constant to match the 
true profile outside the core (at r = 0.2 R Q ). The specific energy profile is computed using 
Eq. |4j Finally, the acceleration of a SPH particle is due either to gravity or to gas pressure. 
These two components are computed using the same particle mass for all particles except 
the core and the companion, for which we distinguish between the gravitational and SPH 
masses. The gravitational mass of the core is M c and its SPH mass is set to balance the 
gravitational acceleration of the envelope and prevent the star from collapsing. As for the 
companion, we treat it as an N-body particle so its SPH mass is M & . 

For the Enzo simulations, the grid is initialized using the stellar model of the primary 
with the addition of a PM particle that represents its core. We fill the computational domain 
with a constant background density to prevent the star from expanding and set the ratio 
between the background density and the minimum density of a cell that belongs to the 
primary to 10~ 4 . This setup is not initially numerically stable. The star tends to expand, 
so we let the initial configuration evolve for a few dynamical times in the absence of the 
companion, while damping the velocity field by a factor of 2 after each cycle. Finally, we 
evolve this relaxed model normally for another few dynamical times to obtain a numerically 
stable model. As a side effect of the relaxation to hydrostatic equilibrium, the Enzo models 
are a little bit bigger — the lower the resolution, the larger the radius of the primary is - 
thus the initial orbital separations between the models are slightly different (Table [I]). 



4. Results 

In this section, we describe the results obtained from our 15 simulations. Since the 
qualitative behavior is the same in all of them, we detail the 0.6 M & case (SPH2, Enzo2 and 
Enzo7). 



4.1. Description of the rapid infall phase 

As explained in £J3j the companion is placed at the surface of the primary. Thus, the 
primary extends beyond its Roche lobe and unstable mass transfer starts immediately. The 
companion, surrounded by stellar matter, exchanges momentum and energy with this gas 
through drag. The orbital separation shrinks on a dynamical timescale and its evolution 
for the 256 3 Enzo simulations is shown in Fig. |4J Although the orbit is initially circular, 
it quickly develops eccentricity due to the geometry of the gas ejection. In order to define 
quantitatively the end of the rapid infall phase and the final orbital separation ad hoc, we 
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consider the evolution of the orbital decay (Fig. pi). As expected, the orbital decay is initially 
quite high (~ 0.01 day -1 ), decreases as less gas is available for the companion to exchange 
energy with, and eventually reaches a plateau. We decide to define the end of the rapid 
infall phase to be at the start of this plateau, which occurs at about 280 days for the 0.6 M Q 
companion (Fig. [5]). All the simulations show the same trend and the lighter the companion, 
the deeper it falls and the longer it needs to reach its final orbital separation. The duration 
of the rapid infall phase is 260, 280, 280, 300 and 340 days, for the 0.9, 0.6, 0.3, 0.15 and 
0.1 M Q companion, respectively. 
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Fig. 4. — Separation between the primary core and the companion as a function of time for 
the 256 3 Enzo simulations. The companion masses are 0.9 (blue), 0.6 (green), 0.3 (red), 0.15 
(cyan) and 0.1 (purple) M Q . 



As orbital energy is transferred to the envelope, the latter is ejected, initially in the 
orbital plane; at later phases there is an almost equal distribution of matter into the polar 
direction as well (Fig. [6]). Overall, almost 90% of the envelope is ejected within an angle of 
30° on each side of the equatorial plane. We compare the orbital velocity of the companion 
(Fig. [7]) with the local sound speed of the gas (Fig. [3j bottom left panel). The former does 
not exceed 50 kms" 1 while the highest sound speed encountered is about 60 kms™ 1 . The 
companion moves only slightly above or below the local sound speed. We therefore conclude 
that the SPH noise could not significantly influence the solution. Also, since the motion of 
the companion is not highly supersonic, the shocks are not strong and we can use Enzo with 
the faster Zeus solver. 
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Fig. 5. — Evolution of the separation (top) and of the orbital decay (bottom) for Enzo7. 
The orbital decay is computed using orbital separations averaged over each cycle (red dashed 
line). The blue vertical line shows the time when we define the end of the rapid infall phase. 



Unlike the SPH computational domain, the Enzo grid is spatially limited. Thus, the 
evolution of the gas that leaves the grid cannot be followed. Therefore, we use the SPH2 
simulation to study the global evolution of the angular momentum and the energy of the 
system. 

We compute the angular momentum using the center of mass of the SPH particles as 
the center of reference. As shown in Fig. [8] for the 0.6 M companion case, the total angular 
momentum of the system is conserved to less than 1%. Since the ejection of the gas is 
asymmetric, the center of reference is eventually located outside the orbit. Consequently, 
studying the orbital components individually is irrelevant as the sign of each component 
changes during a single orbit. Therefore, we study their sum J or b instead. During the 
first 50 days, angular momentum from the orbit almost equally spins up the envelope and 
unbinds mass from the outer layers. Later on, no more additional mass gets unbound (see 



5.1.2) and the angular momentum lost from the orbit spins up the bound envelope only. 



Since the unbound mass is located at large distances from the primary's core, there is no 
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more exchange of angular momentum between the unbound mass and the rest of the system. 
After ~ 150 days, there is no more angular momentum exchange in the system. The primary's 
core and the companion — which are the main contributors to the calculation of the centre of 
mass — switch positions twice per orbit, which leads to small periodic motions of the center 
of mass. These periodic displacements are the causes for the small angular momentum 
fluctuations of the orbital components and the bound mass occurring after 100 days. 

We plot the various energy components in Fig. [9j We start by explaining the different 
components of potential, thermal and kinetic energies represent and how they are computed. 
Among numerous other attributes, each particle i possesses a specific gravitational potential 
energy 0j, a specific thermal energy Ui and a specific macroscopic kinetic energy fcj. By 
definition, 

&= E ~ G ^— (15) 

j particles, j^i ^ 

where G is the gravitational constant, M grav is the gravitational mass of particle j and r^ 
is the distance between particles i and j. We compute these different components using the 
gravitational mass of the particle for the gravitational potential energy and the macroscopic 
kinetic energy, and the SPH mass for the thermal energy (see q3J): 



% = Mr y (/>i (i6) 

Ki = Mf w k % (17) 

Ut = M^ Ul (18) 

where M 4 sp is the SPH mass of particle i used to compute its acceleration due to pressure. 
We recall that both masses are identical for all particles except the primary's core and the 
secondary. Thus, the total gravitational potential energy of the system is 

$tot = 2 E $ < ( 19 ) 

i particles 

Finally, we subtract the contribution of the secondary from the total potential energy 
in order to calculate the binding energy of the envelope: 

$e„v = $tot " Mf a > 2 (20) 
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where the subscript "2" stands for the secondary. 

During the first 200 days when most of the in-spiral happens, the total internal energy of 
the system decreases by more than a factor of two: the envelope expands and therefore cools. 
The energy released is transferred mostly into macroscopic kinetic energy of the gas: the 
envelope is lifted up, accelerated and the outermost part of the envelope becomes unbound in 
the first 50 days. At later times, more energy is transferred from the orbit to the envelope but 
no more material becomes unbound. One can easily note in Fig. [9] how the variations of the 
orbital energy of the core-secondary system and of the total energy of the envelope balance 
each other. The total energy of the envelope remains negative throughout the simulation. 
We follow the evolution of the unbound particles and determine their initial position in the 



envelope. Fig. [10| shows the cumulative mass of the particles that will eventually get unbound 
as a function of their initial distance from the core. It confirms that the unbound mass was 
initially located in the outer part of the envelope and that almost all gas located initially 
closer than 40 R & from the primary's core remains bound at the end of the simulations. 



4.2. Code comparison 

The fact that a code solves the equations in an accurate and precise way in a particular 
situation does not necessary mean it will do so in another regime. Thus, a direct comparison 
of simulations of the CE interaction using two different numerical methods is a good solution 



for testing the ability of the two methods to model this problem. One can see in Fig. 11 and 
Table [T] that for each binary system, the final separations in the Enzo simulations are very 
close to those obtained with the equivalent SNSPH simulations. We may then compare the 
mass evolution of the material in the volume defined by the Enzo grid, the matter within the 
initial volume of the primary, and within the current separation. For the 0.6 M & companion 



(Fig. 12), both the mass within the Enzo grid and the mass within the initial volume of the 
progenitor agree well between the Enzo and the SNSPH runs. For the mass within the orbit 
we notice a difference of ~ 10 -2 M between the Enzo and the SNSPH runs. This difference 
is large compared with the mass of a SPH particle (~ 10 -6 M ) and is due to how accurately 
accretion of the gas by the core and the companion is resolved by the two codes. We have 



plotted, in Fig. 13, density profiles at different times along the line joining the primary core 
and the secondary, for the three simulations with the 0.9 M companion. Accretion onto 
the secondary is better resolved in the SPH simulations in which the maximum density of 
the matter accreted by the companion is about 10 -3 gem -3 . This maximum value depends 
on the resolution of the runs. In the single-grid Enzo runs, accretion is poorly resolved due 
to the low number of cells resolving the local region around each particle. Although mass 



is still accreted around the particles, it eventually becomes dispersed. For the SPH runs, 
around 60 particles interact within a smoothing length so the accretion zone is well resolved. 
On the other hand, the cell width of the Enzo 256 3 runs is about 1.6 R Q so the accretion 
zone cannot be resolved although it is still better than for the Enzo 128 3 simulations as can 



be seen from comparing the different density profiles at 50 days (Fig. 13). However, the 
accurate simulation of accretion onto the secondary is not crucial for the global evolution of 
the system: as we mentioned earlier, the evolution is not driven by accretion but by drag 
forces. Although the density of the matter accreted by the companion differs by up to 3 
orders of magnitudes between the two methods, the accreted mass is negligible compared 
with the companion mass and the final orbital separations are very similar. 



Ricker fc Taam (2008) used the FLASH code (Fryxell et al. 2000) to study the CE 



evolution of a binary system consisting of a 1.05 M & RGB star having a 0.36 M & core and a 
0.6 M & companion. Their implementation is somewhat different from ours since they treat 
the red giant core and the companion as spherical clouds of particles. In spite of those 
differences, their progenitor is almost identical to ours and they find a final separation of 20 
Rq which falls within the range of the results given by our simulations SPH2, Enzo2 and 
Enzo7. Moreover, one can see in Fig. [7] that for the 0.6 M & companion, the velocity of the 
companion stays below 50 kms -1 and therefore, the gas flows are subsonic except in the 



outer layers. This conclusion was also reached by Ricker & Taam (2008) 



4.3. The impact of initial conditions 

In order to determine the sensitivity of the final state of the system to the initial param- 
eters, we start with the Enzo3 simulation and increase by 5% either the initial velocity of 
the secondary (Enzoll) or the initial separation between the two particles (Enzol2), which 
correspond to initial eccentricities of 0.10 (Enzoll) and 0.05 (Enzol2). The evolution of the 



separation for those three simulations is compared in Fig. 14 For Enzoll and Enzol2, the 



ratio of the initial velocity of the companion to the velocity required for a circular orbit is 
higher than one (vo/v C i rc > 1), so the separation must first increase. The larger the orbital 
separation, the more delayed the rapid infall phase is and the later the system reaches its 
final separation. The final separations for Enzo3, Enzoll and Enzol2 are 11.7, 12.0 and 
12.2 R Q , respectively, and the final eccentricities are 0.09, 0.17 and 0.18, respectively. As 
expected, the companion that moves outwards the farthest initially, sinks into the envelope 
with a higher orbital decay velocity. Therefore, it attains a more eccentric orbit and com- 



pletes fewer revolutions around the primary core (Fig. 14). However, the standard deviation 
of the final separation between the three simulations (a ~ 0.2 R Q ) is more than 10 times 
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smaller than the width of a cell. Consequently, we conclude that the final results are quite 
insensitive to the initial conditions at the level tested. 



4.4. Gravitational vs Hydrodynamic drag 

The drag exerted on the companion has two components: gravitational and hydrody- 
namical. The former is due to gravitational forces from matter flowing past the companion 



and colliding with its wake (Bondi & Hoyle 1944 Iben & Livio 1993), while the latter is 



due to ram pressure forces on the companion, 
estimated as: 



The hydrodynamical contribution can be 



^hydro ~ P^l X T\R\ 



(21) 



where R 2 is the radius of the secondary, V2 is the relative velocity between the secondary 
and the envelope, and we have taken the coefficient of drag to be unity for simplicity. In a 



similar manner, the gravitational drag is approximated by (Iben & Livio 1993): 



-Fgrav ~ pv\ X TxR 2 A 



(22) 



where the accretion radius Ra is defined as: 



R- 



2GMo 



(23) 



where c s is the sound speed of the medium. Choosing |v| = 2c s = 80 kins" 1 with an 0.6 M Q 
companion yields Ra ~ 30 R Q . Assuming R 2 ~ 1 -R©, we conclude that the hydrodynamical 
drag is of the order of almost 1 000 times smaller than the gravitational drag, thus negligible. 

This conclusion is also confirmed by the outcomes of our simulations. Indeed, the 
primary's core and the companion are treated as point masses and are not pressure sources, 
except for the primary's core in the SNSPH simulations. Instead of being caused by the finite 
size of the particles, hydrodynamical drag in the models is thus due to the matter accreted 
around them. We pointed out earlier that the accuracy with which accretion was treated 
was different between the two different models because of the different finest resolutions and 
softenings used: accretion is poorly modeled in the Enzo simulations whereas in the SNSPH 
simulations, the companion builds up a sphere of accreted matter about a few R & wide 



around itself (Fig. 13). This should lead to differences in the magnitude of hydrodynamic 
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drag forces. Nevertheless, the consistency of the results suggests that the hydrodynamic 
drag is unimportant in the evolution of the system, confirming the results of |Ricker &: Taam 



(2008). 



Discussion 



5.1. Comparison of simulations and observations 

We now compare the numerical results with a sample of 61 observed post-CE systems 



listed in Zorotovic et al. (2010) and De Marco et al. (2011). 



5.1.1. Final separations 

For a given companion mass (or alternatively mass ratio q) we obtain 3 values for the 
final separation Af, one for each simulation carried out with that companion mass (Table II] 



and Fig. 15 ). One can distinguish between these values at high q (q > 0.34), which correspond 
to "heavy" companions (M 2 > 0.3 M ), and the ones at low q (q < 0.34) corresponding to 
"light" companions (M 2 < 0.3 M ). At high q, the values of Af are very similar and the 
standard deviation is more than 20 times smaller than the average value of Af. At low 
q, the companion sinks deeper and as a consequence, the resolution used in the 128 3 Enzo 
simulations is not sufficient. However, as one increases the resolution to 256 3 cells, the final 
separations converge to the solutions given by the SNSPH simulations. 



Fig. [16] shows the distribution of orbital separations reached by the 61 post-CE systems. 
For all these systems, there has been no substantial orbital shrinkage due to phenomena 



such as magnetic braking or radiation of gravitational waves (see discussion in Schreiber & 



Gansicke 2003). Although they cover a significant range in secondary masses, going from a 



1.1 M & MS star down to a 0.05 M brown dwarf, all of them have separations smaller than 
11 Rq. Furthermore, 87% of those systems have separations smaller than 4 R Q , which is 
smaller than any value obtained in our simulations. This is even more obviously shown in 
Fig. [TrJ, where the final separations for simulations presented here and in the literature are 
compared to the orbital separations of the observed post-CE systems. Although a couple 
of observed systems have q > 0.5, one clearly sees that the simulations with M 2 = 0.9 and 
0.6 M leave the companion far out. Systems with lower mass companions (M 2 < 0.3 M ) 
have by and large lower orbital separations than in our simulations. Simulations of |Sandquist 



et al. (1998) and Ricker & Taam (2008) shown in Fig. 17 give results consistent with ours 



All these numerical simulations suggest that the separations between the secondary and the 
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primary's remnant at the end of the simulated rapid infall phase are too large to explain 
the orbital separation of the currently observed post-CE systems. This suggests that further 
evolution of the orbital separation must occur during the phase immediately following the 



rapid infall phase. We discuss this point further in £5.2 



5.1.2. The state of the envelope at the end of the simulations 

As shown on Table [2j most of the primary's envelope remains bound in all of our 
simulations. We study the situation in detail for our canonical model with the 0.6 M 



companion here. The evolution of the mass for different components is plotted in Fig. 18 It 
first confirms that some envelope mass is unbound only during the first 50 days, after which 
neither angular momentum (Fig. [8]) nor kinetic energy (Fig. [9]) are exchanged between the 
unbound mass and the rest of the system. It also shows that more than 85 % of the mass 
remains bound at the end of the simulation. This outcome, already pointed out by |Sandquist 



et al. (1998), is quite intriguing, since the post-CE binaries observed must have succeeded 



in ejecting their envelope. After about 400 days, most of the envelope mass in our models 



has been moved to a larger radius (~ 100 R Q , see bottom panel in Fig. 19), well outside the 
orbit of the primary core and the companion but remaining bound. 



M 2 (M ) M bound (M G 



SPH1 


0.9 


0.44 


SPH2 


0.6 


0.44 


SPH3 


0.3 


0.45 


SPH4 


0.15 


0.46 


SPH5 


0.1 


0.48 



Table 2: Amount of the envelope mass still bound at the end of the SNSPH simulations. 



a At the start of the simulations, Mb OU nd equals the total envelope mass M e = Mi — M c = 0.49 M 

We now investigate how bound the final system is. We consider the center of mass 
of the system composed by the secondary and the mass within the current orbit as the 
center of our frame of reference. Then, we partition the domain into concentric shells with 
identical thickness, calculate the average radial velocity of each shell and compare it to the 



escape velocity at that location. Fig. [19] shows the escape velocity and the average radial 
velocity of the shells. The radial velocity is always positive and is similar to the space 
velocity at radii larger than 600 R & , as expected for envelope ejection. At radii smaller than 
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300 Rq, the radial velocity is much smaller than the space velocity, suggesting that orbital 
motions dominate at those radii. All the mass within 10 3 Rq is bound, which corresponds 
to more than 85% of the envelope mass. The remaining mass is found at radii between 10 3 
and 6 x 1O 3 _R , where the radial velocities are typically between 25 and 75 kms" 1 . Those 
particles were initially in the outer parts of the giant star, and were the first to encounter the 
secondary. At that time of the in-spiral, the shock was slightly supersonic (V2 ~ 35 kms -1 
and c s ~ 20 kms -1 ). This regime of evolution is thus different from later phases when 
the secondary sinks deeper into the primary's envelope, where its velocity does not really 
increase (Fig. Pf]) but the sound speed of the medium does (Fig. [3]). 



We can measure how much extra energy would be required to unbind the envelope at 
each radius, using the definition 

^extra = ^ ~M t (v e>i - V r rf (24) 

i 

where v e ^ and v r ^ are the escape velocity at the location of the i-th shell and its average 
radial velocity, respectively. One finds -E ex t ra ~ 8.4 x 10 45 ergs which represents just over 
10% of the initial binding energy of the primary envelope. Thus, a relatively small additional 
input of energy could be sufficient to completely unbind the remaining envelope material. 

We have compared here the final separations deduced from observations and those de- 
termined from the simulations. We have purposefully stayed away from calculating the 



ejection efficiency a (Webbink 1984 De Marco et al. 2011). Indeed, we question what the 



relevance of calculating a is when the envelope has not yet been fully ejected, true both in 



the Sandquist et al. (1998) and our simulations. We therefore defer for the moment the task 
of calculating a from simulation — a long term goal of this project — until the simulations 
are more advanced. 

In conclusion, the hydrodynamic simulations do not reproduce the post-CE systems in 
the sense that the system is left at too large separations and the envelope is not unbound 
at the end of the rapid infall phase. This means that either physical processes that are not 
accounted for in the simulations are responsible for the envelope ejection, or the envelope 
ejection and a significant reduction of the orbit actually happens during the later subsequent 
slow in-spiral phase. We discuss both possibilities in the following section. 
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5.2. Reproducing the observations 

In this section we first study and quantify physical processes that are not taken into 
account in our hydrodynamic simulations and that might be responsible for ejecting the 
envelope. Then, we focus on the subsequent slow in-spiral phase and investigate whether 
the envelope can be ejected and the separation significantly reduced during this subsequent 
phase. 



5.2.1. Rotation of the primary 

The envelope of the progenitor is initially non-rotating and although the calculation done 
in §|3] shows that, regardless of the initial rotation velocity of the envelope, its rotational 
energy is negligible in comparison with its binding energy, we suspected at first that the 
absence of rotation might be the reason for most of the envelope to remain bound. However, 



Sandquist et al. (1998) carried out two identical simulations where they modified the initial 
rotation state of the primary from a giant star in synchronization with the orbit to a non- 
rotating one (their simulations 1 and 2). In both cases, the evolution of the bound mass 
and the final orbital parameters are similar. It thus does not seem that changing the initial 
rotation of the primary leads to a different CE outcome. 



5.2.2. Physics not included in the simulations 
The hydrodynamics codes use an ideal gas equation of state ( §2TJ ) which, by definition, 



does not include variable abundances and the different ionization layers of the envelope. Han 



et al. ( 1995 ) suggested that recombination might play a role in CE interactions. As the outer 
parts of the envelope expand and cool, ions recombine with electrons, releasing energy that 
could aid in unbinding the envelope. Although it is unclear how efficient this process is and 
how much of the initial recombination budget can be used, one can calculate an upper limit 
on how much energy can be injected into the envelope by recombination. 

According to our stellar evolution model, the hydrogen fraction within the convective 
envelope of our RGB star is X ~ 0.68. The mass of the envelope is M e = 0.49 M & and 
each proton recombining with an electron produces an energy E = 13.6 eV. We also have to 
calculate how much of the envelope is ionized. Therefore, we calculate the partition functions 
Z for hydrogen. The hydrogen ion has no degeneracy so Z<i — 1. The partition function for 
the hydrogen atom at temperature T is 
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Zi = y^ 2n 



2 E (l/n 2 - 1) 
exp - 



n=l 



A;T 



(25) 



where k = 8.6173 x 10~ 5 eV.K^ 1 is the Boltzmann constant. We truncate the sum in Eq. 25 
at the first integer n max such that the distance at which the electron orbits the proton for this 
quantum number is larger than Z max = 10~ 6 cm, i.e. a n^ ax > / max where a = 5.2918 x 10 9 cm 



is the Bohr radius (Miranda 2001). We then use the Saha formula to calculate the ratio of 



ionized to neutral hydrogen (Carroll & Ostlie 1996): 



*/*=ap£T-<-^ 



(26) 



where n e is the number density of free electrons and m e is the electron mass. We find 
that 91% of the envelope is ionized. Consequently, the recombination of the whole ionized 
envelope would produce an extra energy 



R 



recomb 



0.91 x IE 



N t 



M. 



x 13.6 eV 



(27) 



H 



where Na is the Avogadro number and Mh the atomic mass of hydrogen. One finds -Erecomb ~ 
1.18 x 10 46 ergs, which is slightly higher than the extra energy E extTa required to eject the 
envelope in our canonical model (£5.1.2). Thus, we conclude that recombination in the 



envelope could substantially aid in unbinding it. 

Another source of energy could be radiation pressure. For low- and intermediate-mass 
giants in hydrostatic equilibrium, radiation pressure (P ra d = aT 4 /3, where a is the radiation 
constant) is negligible compared to gas pressure (Eq. lib: for our primary, P ra d/ Pgas i$ 0.01 
except in a small zone (0.1 -R Q < r < 10 R Q ), where P ra dl Pgas % 0.1. However, the deep 
in-spiral of the companion within the primary's envelope will induce local shock heating. 



The increase of temperature is proportional to the square of the Mach number (Tarbell et al. 



1999), so even if the companion is orbiting at twice the local sound speed, the radiation 



pressure to gas pressure after the shock becomes: 



P 



rad 



gas 



after 



OC 



P 



rad 



P 



gas 



before 



(M : 



2\3 



6.4 



(2* 



Therefore, including radiation pressure in the equation of state will increase the total pressure 
locally and might reduce the energy required to eject the envelope. However, it is possible 
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that this effect is globally small, since this extra heating source is probably very localized 
around the companion. 



5.2.3. The post-rapid-infall phase 

At the end of the rapid infall phase, the orbit is stable until the end of the simulations 
(a few more years). Consequently, there is no further hydrodynamical coupling between the 
extended envelope and the surviving binary. We now investigate whether the envelope is 
likely to be ejected during this slower in-spiral phase. 

Although the resolution of the simulations prevents us from quantifying how much 
envelope will be left around the core of the primary, one can still describe qualitatively what 



the evolution of the primary's remnant will be. Fig. 18 shows that less than 10 2 M is left 



around the primary's core, so the primary will depart the giant branch (Bloecker 1995, but 



see also the discussion in De Marco et al. 2011). Then two scenarios might occur depending 



on how long the partially ejected envelope will take to fall back. 

If the star is given enough time to transit to the blue due to hydrogen burning at the 
base of the envelope before the lifted envelope falls back, the star will readjust on its thermal 
timescale of the remaining envelope, and eventually end its life as a Helium white dwarf. 
This transition will last ~ 10 3 years during which the star will have a luminosity between 



300 and 1000 L e (Iben & Tutukov 1993, their Fig. 1), which is consistent with the more 
recent work of Driebe et al. ( 1998[ ) (their Fig. 1). If we assume the remnant to have a 
luminosity L c ~ 500 L Q , we can compare the gravitational acceleration of a gas particle 
with the radiation acceleration defined by 



H 



^rad 



Anr 2 c 



(29) 



where r is the distance between the gas particle and the core, k = 0.4 cm 2 g _1 is the opacity 
for Thompson scattering for hydrogen, and c is the speed of light. We still find the radiation 
acceleration to be overall almost two orders of magnitude smaller than the gravitational 
acceleration. 

If on the contrary, the envelope falls back before the primary's remnant had crossed the 



Hertzsprung-Russell diagram, a circumbinary disk will form (Kashi & Soker 2011). They 



refer to the numerical work done by Artymowicz et al. (1991), which suggests that in such a 
configuration, the binary separation will decrease due to Lindblad resonances — mainly - 
as well as viscous tides. Although this mechanism has the advantage of explaining how the 
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orbital separation will diminish during the subsequent phase, the ability of radiation to eject 
the gas will even be reduced in comparison with the previous situation, so it is not clear how 
the latter will eventually be unbound. 

In conclusion, radiation acceleration alone does not seem to be responsible for unbinding 
the remaining gas, regardless of the time the partially ejected envelope will remain suspended 
for. 



6. Summary and Future Work 

In this work we have carried out three-dimensional hydrodynamic simulations of the 
CE interaction between a 0.88 M Q RGB star and companions with mass ranging from 0.1 
to 0.9 Mq. We have used both an Eulerian grid code (Enzo) and a Lagrangian SPH code 
(SNSPH) with various resolutions. They both have advantages and disadvantages and can 
be used for different purposes: while one might rather use SPH to study the accretion 
around the secondary, even a uniform grid code is more suitable in resolving the low-density 
extended envelope. Of course, adaptive mesh refinement combines the advantages of both 
of these methods at the cost of increased code complexity. 

We first compared the outcomes of those simulations with each other. We found that 
the results are very similar for companion masses Mi > 0.3 M . We thus conclude that in 
this regime, the resolutions used are sufficient to study the global evolution of the system 
during the rapid infall phase of the interaction, which is driven mainly by gravitational drag. 
For lower companion masses (M 2 < 0.3 M ) that penetrate deeper in the giant's envelope, 
the 128 3 Enzo runs are under-resolved but the Enzo results converge to the solutions from 
the SNSPH simulations. 

We then compared the outcomes of our simulations with observed post-CE systems. The 
final separations are found to be systematically higher than those deduced from observations, 



as is the case for the past simulations by Sandquist et al. (1998), De Marco et al. (2003) 



and Ricker & Taam (2008). Moreover, mass is only unbound during the early stages of 
the interaction (~ 50 days for the 0.6 M companion) and most of the envelope remains 
bound at the end of the simulations, as was the case for the earlier simulations of |Sandquist 



et al. (1998). We investigated whether there might be additional processes that were not 
accounted for in the simulations. We found that recombination can contribute significantly, 
but stellar rotation and radiation pressure play only marginal roles. Finally, we wondered 
whether the bound envelope is a result of imprecise simulations or a real physical feature. If 
the latter, then one would have to follow the subsequent evolution of the system to determine 
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the actual outcome of the CE. Fall back disks may form and even have an impact on the 



inner binary (Artymowicz et al. 1991, Kashi & Soker 2011) 



After the submission of this paper, Ricker and Taam made their paper Ricker fc Taam 



(2011) available. This paper continues the work introduced in Ricker & Taam (2008). In 



their simulation, only about 25 % of the primary's envelope is unbound. Although this value 
is slightly higher than ours, it is in agreement with our work in the sense that most of the 
envelope remains bound. They also claim that the ejection occurs mostly in the orbital 
plane, as it is the case in our simulations. However, the extended envelope at the end of 
their simulation is rotating much faster than it is expanding which is in contradiction with 
our results 



5.1.2) but might be due to the fact that their primary is initially rotating. 
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Fig. 6. — Density slices in the orbital plane (left) and in the perpendicular plane (right) at 
0, 50, 85 and 130 days (from top to bottom) for the Enzo7 simulation. The scale used for 
the velocity vector field is the same on each frame and is such that the velocity shown on 
the top panel equals the initial orbital velocity of the primary (~ 23 kms -1 ). 
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Fig. 7. — Evolution of the companion velocity for the Enzo7 simulation. 
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Fig. 8. — Evolution of the z-component of the total angular momentum (Jtot), the angular 
momentum of the core and the companion (J or b), the angular momentum of the bound mass 
(Abound) and the angular momentum of the unbound mass (./unbound) f° r the SPH2 simulation. 
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Fig. 9. — Energy components for the SPH2 simulations. Plotted are the total energy (E tot ), 
the total gravitational potential energy $tot, the internal energy of the system (Utot), the 
gravitational potential energy of the envelope ($ env ), the gravitational potential energy from 
the core-companion interaction ($ C 2), the kinetic energy of the core and the companion 
(K c2 ), the kinetic energy of the bound mass (K h ), the kinetic energy of the unbound mass 
(K n ), the orbital energy of the core-companion system (E c2 ) and the total energy of the 
envelope (E env = $ cnv + U to t + K\>)- The beat frequency seen on K c2 and <J> C 2 are due to the 
non-synchronization between the orbital period and the data dumping frequency. 
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Fig. 10. — Initial distribution within the envelope of the mass that will eventually get 
unbound for SPH2. 
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Fig. 11. — Separation between the core of the primary and the 0.6 M companion as a 
function of time for the SPH2 (left), Enzo2 (middle) and Enzo7 (right) simulations. Again, 
the beat frequency seen in the SPH simulation is due to the non-synchronization between 
the orbital period and the dumping frequency. 
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Fig. 12. — Each panel shows the mass within the equivalent Enzo grid (plain), the inital 
volume of the primary (dash) and the orbit (cross-solid) as a function of time for the SPH2 
(left), Enzo2 (middle) and Enzo7 (right) simulations. 
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Fig. 13. — Density profiles along the line joining the core and the 0.6 M companion at 
(dotted line), 50 (dash-cross line), 100 (dash-dot line), 300 (dashed line) and 500 (solid 
line) days for SPH2 (top), Enzo2 (middle) and Enzo7 (bottom). The vertical lines show the 
position of the companion. 
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Fig. 14. — Left: separation between the core of the primary and the companion as a func- 
tion of time for the Enzo3 (solid blue), Enzoll (dashed red) and Enzol2 (dash-dot green) 
simulations. Right: a detail of the comparison from the left panel at ~ 340 days. 
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Fig. 15. — Final separations as a function of the mass ratio q for the SNSPH (black cross) 
Enzo 128 3 (blue circle) and Enzo 256 3 (red triangle) simulations. 
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Fig. 16. — Distribution of post-CE systems as a function of their observed orbital separation 



from Zorotovic et al. (2010) and De Marco et al. (2011). 
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Fig. 17. — Comparison between the orbital separations of observed post-CE systems (black 
dot) and the final separations reached at the end of the simulations (red circle), as well as the 
ones by Sandquist et al. ( 1998[ ) (green circles) and by Ricker k, Taam| ( 2008 ) (blue triangle). 
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Fig. 18. — Evolution of the total mass (solid line), the bound mass (dashed line), the mass 
within the volume of the Enzo grid (dotted line), the mass within the initial volume of the 
primary (dash-cross line) and the mass within the orbital separation (dash-dot line) for the 
SPH2 simulation. 
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Fig. 19. — Top: Comparison between the escape velocity (dashed line) and the radial velocity 
(black dots) of the final system for the SPH2 simulation. Bottom: Mass enclosed as a function 
of radius. 



